Relic gravitational waves from light primordial black holes 



Alexander D. Dolgov a ' &,c , Damian Ejlli a ' 6 
January 20, 2013 

y—i ■ 

^ ■ a Universita degli Studi di Ferrara, 1-44100 Ferrara, Italy 

■ b INFN, Sezione di Ferrara, 1-44100 Ferrara, Italy 

3 ! C ITEP, 113259 Moscow, Russia 

i— i ■ 

Abstract 

The energy density of relic gravitational waves (GWs) emitted by primordial black holes (PBHs) 
is calculated. We estimate the intensity of GWs produced at quantum and classical scattering of 
PBHs, the classical graviton emission from the PBH binaries in the early Universe, and the graviton 
emission due to PBH evaporation. If nonrelativistic PBHs dominated the cosmological energy 
density prior to their evaporation, the probability of formation of dense clusters of PBHs and their 
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binaries in such clusters would be significant and the energy density of the generated gravitational 



waves in the present day universe could exceed that produced by other known mechanisms. The 
intensity of these gravitational waves would be maximal in the GHz frequency band of the spectrum 
or higher and makes their observation very difficult by present detectors but also gives a rather good 
possibility to investigate it by present and future high frequency gravitational waves electromagnetic 
j> \ detectors. However, the low frequency part of the spectrum in the range / ~ 0.1 — 10 Hz may be 

C*~) ■ detectable by the planned space interferometers DECIGO/BBO. 

For sufficiently long duration of the PBH matter dominated stage the cosmological energy 



{Sj \ fraction of GWs from inflation would be noticeably diluted. 
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1 Introduction 



Since the prediction of gravitational waves (GW) by Albert Einstein in 1916 [1J on the basis of general 
relativity, they have been an object of intensive studies. Gravitational waves are thought to be fluc- 
tuations in the curvature of space-time, which propagate as waves, traveling outward from the source. 
Although gravitational radiation has not yet been directly detected, it has been indirectly shown to exist 
because it increases the pulsar orbital frequency [2] in good agreement with theoretical predictions. 

Roughly speaking there are two groups of possible sources of gravitational radiation which may be 
registered by gravitational wave detectors either on the Earth or by space missions. The first group 
includes energetic phenomena in the contemporary universe, such as emission of GWs by black hole 
or compact star binaries, supernova explosions, and possibly some other catastrophic phenomena. The 
second group contains gravitational radiation coming from the early Universe, which creates today an 
isotropic background usually with rather low frequency. Such gravitational radiation could be produced 
at inflation, phase transitions in the primeval plasma, by the decay or interaction of topological defects, 
e.g. cosmic strings, etc. 

The graviton (gravitational wave) production in the Friedmann-Robertson- Walker metric was first 
considered by Grishchuk [3J, who noticed that the graviton wave equation is not conformal invariant 
and thus such quanta can be produced by conformal flat external gravitational field. Generation of 
gravitational waves at the De Sitter (inflationary) stage was studied by Starobinsky [3] (see also ref. [5]). 
The stochastic homogeneous background of the low frequency gravitational waves is now one of the very 
important predictions of inflationary cosmology, which may present a final proof of inflation. 

In this work we discuss one more source of gravitational wave (GW) radiation in the early Universe, 
namely, the interaction between primordial black holes (PBH). We consider relatively light PBH, such 
that they evaporated before the big bang nucleosynthesis (BBN) and so they are not constrained by 
the light element abundances. Cosmological scenario with early formed and evaporated primordial 
black holes producing gravitons was considered in ref. [5J. Here we will remain in essentially the same 
frameworks and study in addition the GW emission in different processes with PBH. 

According to ref. [8] the life-time of evaporating black hole with initial mass M is equal to: 

10240 7T M 3 . . 

TBH = -T7 -T- , (1) 

J\ e f f m Pl 

where the Planck mass is mpi = 2.176 ■ 10~ 5 g and N e ff is the number of particle species with masses 
smaller than the black hole temperature: 

T BH = ™*L . (2) 

To avoid a conflict with BBN the black holes should had been evaporated before cosmological time 
t ~ 10~ 2 s [9] and thus their mass would be bounded from above by 

M < 1.75 • 10 8 g. (3) 

The temperature of such PBHs should be higher than 3 • 10 4 GeV and correspondingly N e ff > 10 2 . On 
the other hand, as is discussed in what follows, the PBH mass is bounded from below e.g. by equation 
(16). This is the mass range of PBHs considered in this work. Such PBH are not constrained by any 
astronomical data, which are applicable to heavier ones [5]-[10j. 
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Primordial black holes should interact in the early Universe creating gravitational radiation. Below 
we estimate the efficiency of GW emission in several processes with PBH. In sec. 2 some mechanisms of 
PBH production and PBH evolution in the early Universe are briefly described. We stress, in particular, 
a very important role played by the clumping of PBH due to gravitational instability at the matter 
dominated stage. In section 3 we consider the initial interaction between the PBHs when they started 
to "feel" each other and accelerate with respect to the background cosmological expansion. In sec. 4 
the quantum bremsstrahlung of gravitons at PBH collisions is discussed, which is quite similar to the 
electromagnetic bremsstrahlung at Coulomb scattering of electrically charged particles. Next, in sec. 
5 we consider the classical emission of GW at accelerated motion of a pair of BHs in their mutual 
gravitational field. In sec. 6 we evaluate the energy loss of PBHs due to their mutual interaction. It 
may be relevant to the estimation of the probability of formation of PBH binaries. The gravitational 
radiation from PBH binaries in high density clusters is discussed in sec. 7, In sec. 8 we calculate 
the present day energy density of gravitons produced at PBH evaporation. In sec. 9 we review some 
mechanisms of the production of stochastic background of GWs. In sec. 10 the status of existing and 
planned detectors of GWs is discussed. In sec. 11 we conclude. 



2 Production and evolution of PBH in the early universe. 

Formation of primordial black holes from the primordial density perturbations in the early Universe was 
first considered by Zeldovich and Novikov [11] and later by Hawking and Carr [121 03] • PBHs would be 
formed when the density contrast, Sp/p, at horizon was of the order of unity or, in other words, when 
the Schwarzschild radius of the perturbation was of the order of the horizon scale. If PBH was formed 
at the radiation dominated stage, when the cosmological energy density was p(t) = 3m 2 pl /(32Trt 2 ), and 
the horizon was lh = 2t, the mass of PBH would be: 

M(t) = m 2 Pl t ~4- 10 38 ( — J g (4) 

where t is the time elapsed since Big Bang. 

The fraction of the cosmological energy density of PBH produced by such mechanism depends upon 
the spectrum of the primordial density perturbations. We denote this fraction fl p and take it as a free 
parameter of the model. The data on the large scale structure of the Universe and on the angular 
fluctuations of the cosmic microwave background radiation (CMB) show that the spectrum of the pri- 
mordial density fluctuations is almost flat Harrison- Zeldovich one. For such spectrum the probability 
of PBH production is quite low and f] p < 1. However, the flatness of the spectrum is verified only for 
astronomically large scales, comparable with the galactic ones. The form of the spectrum for masses 
below 10 10 g is not known. Inflation predicts that the spectrum remains flat for all the scales but there 
exist scenarios with strong deviation from flatness at small scales. In particular, in ref. [T4"| [T5] a model 
of PBH formation has been proposed which leads to log- normal mass spectrum of the produced PBH: 



dN 

—— = C exp 

dM 



(M - M ) 
M 2 



(5) 



where C, Mo, and Mi are some model dependent parameters. Quite naturally the central value of PBH 
mass distribution may be in the desired range Mq < 10 9 g. In this model the value of Q p may be much 
larger than in the conventional model based on the flat spectrum of the primordial fluctuations. We will 
not further speculate on the value of fl p and on the form of the mass spectrum of PBH. In what follows 
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we assume for an order of magnitude estimate that the spectrum is well localized near some fixed mass 
value and that Q p is an arbitrary parameter. Different mechanisms of PBH production are reviewed e.g. 
in ref. HE]. 

We assume that PBHs were produced in radiation dominated (RD) Universe, when the cosmological 
energy density was equal to 

^ m Pi ( R \ 

If we neglect the PBH evaporation and possible coalescence, their number density would remain constant 
in the co-moving volume, n B H{t)a 3 (t) = const. In what follows the instant decay approximation for 
evaporation is used. The cosmological evolution of PBHs with more realistic account of their decay was 
studied in ref. [IB] . 

Since the black holes were non-relativistic at production, their relative contribution to the cosmo- 
logical energy density rose as the cosmological scale factor, a(t) : 

n BH (t) = a p , (7) 

where a p is the value of the scale factor at the PBH production and at RD-stage a(t)/a p = (t/t p ) l l 2 . 
The moment t p of the black hole production is connected with the PBH mass through eq. (4). Hence 

M 

t, 



P 2 
m Pl 



Thus if PBHs lived long enough, they would dominate the cosmological energy density and the Universe 
would become matter dominated at t > t eq , where 

M r g 

and r g = 2M/m 2 Pl is the gravitational (Schwartzschild) radius of a black hole. 

In what follows we assume that all PBHs have the same mass M, but the results can be simply 
generalized by integration over the PBH mass spectrum. 

Evidently at RD stage the number density of PBHs drops as: 

n BH (t) = n p (-%] = n p (^] ' , (10) 



a(t) 



n BH (t) = n p (^Y 2 (t-f) . (11) 



while at MD stage 



, teq J V ' 

Cosmological mass fraction of BH as a function of time behaves as 

^ / x n BH (t)M 167T 9 . , , . 

n BH (t) = K ' = r g t 2 n BH {t) , (12) 

Pc <J 

i.e. VL BB ~ t 1 / 2 at RD stage. After the onset of the PBH dominance, VL BB approached unity and 
remained constant till the PBH evaporation when Q BB quickly dropped down to zero and the universe 
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became dominated by relativistic particles produced by PBH evaporation. All relics from the earlier 
RD stage would be diluted by the redshift factor it eq /T BH ) 2 ^ . In particular the energy density of GWs 
produced at inflation would be diminished by this factor with respect to the standard predictions. Such 
dilution may cause problems with baryogenesis. However, these problems may be resolved if baryogenesis 
took place at the process of PBH evaporation through the mechanism suggested by Zeldovich pj5] and 
quantitatively studied in ref. [201 121] • Somewhat similar model of baryogenesis by heavy particle decay 
(e.g. by bosons of GUT) created at PBH evaporation was considered in ref. [221 1231 [2U I2"5] . 

To survive till equilibration the PBHs should live long enough so that their evaporation time t ev 
would be larger than t eq or tbh > t eq — t p which can be translated into the bound on the PBH mass: 

M > (iL) m „ ('-X' 2 , 5.6 • 10- m) ^ (13) 



v 3.2- 10V V^p / V 100 / fi p 

where Q p <C 1 and M is mass of PBHs at production . Both constraints (3) and (13) would be satisfied 
if 

^>0.7-10^(^y /6 . (14) 

For example, if Q p = 10 -10 , the black holes should be heavier than 1.2 • 10 4 g. 

When the Universe became dominated by non-relativistic PBHs, primordial density perturbations, 
A = Sp/p, should rise as the cosmological scale factor. They could reach unity at cosmological time ti 
satisfying the condition: 

// \ 2/3 

A,„(£) ~1, (15) 

where A in is the initial magnitude of the primordial density perturbations. To be more accurate, the 
evolution of density perturbations depends upon the moment when they cross horizon, see below, eq. 
(19). For the moment we neglect this complication to make some simple estimates. 

The initial density contrast is usually assumed to be of the order of A in ~ 10~ 5 — 10~ 4 which is not 
necessarily true at small scales and may be much larger, especially in the model of ref. [T4"l ll~5] . 

Evidently the BH life-time, tbh, must be long enough, so that the density fluctuations in BH matter 
would rise up to the values of the order of unity. The condition t ev > t\ or equivalently tbh > t± — t p 
leads to the following restriction on the PBH mass: 

We can see that eq. (16) puts a stronger lower limit on PBHs mass than eq. (13). The limits are 
comparable only if A in w 1. Using eqs. (16) and (3) we get a stronger than (14) restriction on fl p : 



1 In fact in equation (13) there must be the PBHs mass at the equilibrium time, M(t eq ). Due to evaporation the PBH 
mass as a function of time is given by M(t) = M(t p )(l — I/tbh) 1 ^ and it is easy to see that for tbh > teq it gives 
M = M(t p ) ~ M(t eq ), so hereafter we refer to M as the mass of PBH at production. 



4 



After A reached unity, the rapid structure formation would take place and high density clusters of 
PBHs would be formed. As we see in what follows, generation of gravitational waves would be especially 
efficient from such high density clusters of primordial black holes. 

Let us assume that the spectrum of perturbations is the flat Harrison- Zeldovich one and that a 
perturbation with some wave length A crossed horizon at moment U n . The mass inside horizon at this 
moment was: 

M b (t in ) = m 2 Pl t in . (18) 

It is the mass of the would-be high density cluster of PBHs. This initial time is supposed to be larger 
than t eq (9), i.e. the horizon crossing took place already at MD-stage. For flat spectrum of perturbations 
density contrast, A = Sp/p, at horizon crossing is the same for all wave lengths. After horizon crossing 
the perturbations would continue to grow up as the scale factor, A(t) = A in (//tj n ) 2//3 . Such rise would 
continue till moment ti(tj n ) such that: 

A[ti(t in )] = A m [tx(t m )A m ] 2/3 = 1 or ti(tin) = t in Ar 3/2 . (19) 

The radius of the PBH cluster rose almost as the cosmological scale factor till t = t\(ti n ). After the 
density contrast has reached unity the cluster would decouple from the common cosmological expansion. 
In other words, the cluster stopped to expand together with the universe and, on the opposite, it would 
begin to shrink when gravity takes over the free streaming of PBHs. So the cluster size would drop down 
and both Ubh and p b would rise. The density contrast would quickly rise from unity to A b = Pb/Pc ^ 1, 
where p c and p b are respectively the average cosmological energy density and the density of PBHs in 
the cluster (bunch). It looks reasonable that the density contrast of the evolved cluster could rise up 
to A = 10 5 — 10 6 , as in the contemporary galaxies. After the size of the cluster stabilized, the number 
density of PBH, ubh, as well as their mass density, Pbh, would be constant too. But the density 
contrast, A b would continue to rise as {t/ti) 2 because p c drops down as 1/t 2 . From time t — t\ to 
t = tbh the density contrast would additionally rise by the factor: 

A( TM) = A( tl )(^) 2 =A( t ,)(£) 4 , 

where t\ and Mi ow are given by eqs. (15 ) and (16) respectively. 
The size of the high density clusters of PBH would be 

D _ A -1/3,2/3,1/3 

and the average distance between the PBHs in the bunch can be estimated as: 

* = (M/M b ) 1/3 R b = A- 1/3 /f r 1 / 3 = 2- 2 /3 A - 1 /3 A -ifi p -V3 rg 

It does not depend upon t in . Here eqs. (15) and (9) have been used. 
The virial velocity inside the cluster would be 

So PBHs in the cluster can be moderately relativistic. 
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(20) 



(21) 



(22) 



(23) 



Later, when t = tbh, black holes would decay producing relativistic matter and the Universe would 
return to the normal RD regime. However, the previous history of the earlier RD stage would be 
forgotten. 

For the future discussion it is convenient to introduce the average distance between the PBHs at 
arbitrary time, d = n BI ^ , where ubh = Pbh/M is the number density of PBHs. Since 



n = — 



32ntlMn p 
2>m 2 Pl 



32vr ft 



d. 



the average distance between PBHs at the production moment is equal to 

d p = (4n/3) 1 / 3 r g Q p 1 / 3 . 



(24) 



(25) 



When the mutual gravitational attraction of PBH may be neglected, d rises as cosmological scale factor, 
a(t). 

Gravitational waves produced in the early universe will be hopefully registered in the present epoch. 
The sensitivity of GW detectors strongly depends upon the frequency of the signal. The frequency /* 
of GW produced at time t* during PBH evaporation, is redshifted down to the present day value, /, 
according to: 



/ = /. 



a 



0.34/, 



100 

9s(T* 



1/3 



(26) 



where T = 2.725 K [26] is the temperature of the cosmic microwave background radiation at the present 
time, T* = T(t*) is the plasma temperature at the moment of radiation of the gravitational waves, and 
gs(T*) is the number of species contributing to the entropy of the primeval plasma at temperature T*. 
It is convenient to express T in frequency units, T = 2.7 K = 5.4 • 10 10 Hz. 

The temperature of the primeval plasma after the PBH evaporation can be approximately found 
from: 



P 



m 2 Pl _ n 2 g*(T*)T* 
6rtt 2 " 30 



(27) 



where g*(T*) f=s 10 2 is the contribution of different particle species to the energy density at temperature 
T* and t\ < t < t ev . For relativistic plasma g*(T) = gs(T). Since t ev = tbh + t P — tbh, we obtain from 
equation (27) at time £* = tbh- 



T*(T BH ) : 

Substituting the numbers we find: 

T^Tbh) « 0.01 Imp, 



30 



1/4 



N, 



eff 



6vr 3 ^(T,)J V3-2-10V M 3 / 2 ' 



1/2 5/2 



PI 



(21 



100 

.9s(T*) 



1/4 'n k „y* (u*ly /2 



100 



V M ) 



(29) 



For comparison at the PBH production moment the temperature of the primeval plasma was: 



T w 
± p ~ 



0.2m„ (^) V2 



(30) 
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Using eqs. (26) and (29), we find that the present day frequency of the GWs, emitted at T* (28) with 
frequency /*, would be equal to: 



/ = 1.7- 10 12 Hz 



100 

9 s (T* 



1 1/12 / 100 \ 1/2 / f \ f M \ 3/2 

\ ( h \ I \ (31) 



N eff) \mpi) \m Pl ) 



If we take the maximum frequency of the emitted gravitons f ma x* ~ r g 1 = "^p;/2M, the GW maximum 
frequency today would be: 

/ M \ 1/2 / /If \ 1/2 



« 8.6 • lO^Hz ( — j = 5.8-10">Hz(^- ] (32) 



3 Onset of GW radiation 

Once PBHs enter inside each other cosmological horizon 2 they start to interact and thus to radiate 
gravitational waves due to their mutual acceleration. The corresponding time moment th is determined 
by the condition 2t h = d{t h ) and remembering that it happened still at RD stage, we find 

fc = K^) 2/ W /S - (33) 



2 V 3 



For t > th, the curvature effects can be neglected and the PBH motion is completely determined by the 
Newtonian gravity: 

M BH r 

r = - 34 

m pi r z r 

with the initial conditions = |rj| = d(ti) and |rj| = H(ti)\ri\, where r is the position vector of 
PBHs. For ti = th their relative initial velocity |r-j| = Vi = 1 and non-relativistic approximation is 
invalid. To avoid that we should choose ti > t h such that 1. The solution of the equation of motion 
demonstrates that the effects of mutual attraction at this stage and production of GW are weak. 

After PBHs enter inside each other horizon and Newtonian gravity can be applied, their acceleration 
toward each other becomes essential when their Hubble velocity drops below the capture velocity. The 
corresponding time moment, t c , when it happened, is determined from the condition: 

\aq = \m M ( tc) f < - J^. (35) 

If it took place at the RD regime, the corresponding time moment would be equal to: 

«. = -r?f' (36) 

and the density parameter of PBHs at t — t c would be 

ft\ 1 ' 2 An 

siBH(t e ) = n p [f) =y >1 - (3 7 ) 



2 The cosmological horizon is the distance which PBHs started interacting with each other exchanging gravitons and 
should not be confused with the black hole event horizon. 
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Thus at t — t c the universe is already matter dominated and we have to use the non-relativistic expansion 
law, a ~ t 2 / 3 , starting from the moment t = t eq (9). Accordingly the average distance between BHs, 
when t > t eq , grows as: 

ft \ 1/2 ft \ 2/3 

Now we find that the condition that the Hubble velocity, Vh = (2/3t c )d c is smaller than the virial one, 
for average values, reads: 

Adl 

p < 1 . (39) 



Q ,3/2,1/2 



One can see that this condition is never fulfilled. However, this negative result does not mean that the 
acceleration of BHs and GW emission are suppressed, because of the mentioned above effect of rising 
density perturbations. 



4 Bremsstrahlung of gravitons. 

PBH scattering in the early Universe should be accompanied by the graviton emission almost exactly 
as the scattering of charged particles is accompanied by the emission of photons. The cross-section of 
the graviton bremsstrahlung in particle collisions was calculated in ref. [27] for the case of two spineless 
particles (here black holes) with masses m and M under assumption that m <C M. In non-relativistic 
approximation, p 2 <C m 2 , the differential cross section reads: 



64M 2 m 2 d£ 

da = 7, 

15mJ i 



(40) 



where £ is the ratio of the emitted graviton frequency, u = 2ir f, to the kinetic energy of the incident 
black hole, i.e. £ = 2mu/p 2 . We will use expression (40) for an order of magnitude estimate assuming 
that it is approximately valid for arbitrary m and M, in particular, for m ~ M. 

The energy density of gravitational waves emitted at the time interval t and t + dt in the frequency 
range u and u + du> is given by 

-^ r = Vrein BH U^-)dt, (41) 

where ubh is the number density of PBH and v Te \ is their relative velocity. 

The energy emitted in the frequency interval lo G [0,u; max ] per unit time is proportional to the 
integral 



2 £,max 



5^1^? +-(2-0 In 



3/o ,m 1 + v 7 !^ 



2 V l-VW 



(42) 



The maximum value of the frequency of the emitted gravitons should be smaller than either the 
kinetic energy of the colliding BHs, E^in = p 2 /(2M) or the BH inverse gravitational radius, l/r fl = 
m 2 pl /2M 1 depending on which of the two is smaller. Their ratio is Ek in r g = M 2 v 2 /m 2 pl , so for M < 
mpiv~ x the maximum frequency would be the PBH kinetic energy and in this case £ max = 1. It 
corresponds to the situation when PBH is nearly captured. It looses practically all its kinetic energy, 
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which goes to the graviton. For PBHs in the high density clusters, when v ~ 0.1, the maximum frequency 
would be 0J max ~ l/r g for all PBHs heavier than lOrripi. In this case (, max = (m P i/Mv) 2 . 

The first, rather exotic case, when M < mpi/v can be realized only if Q p > 0.01, see eq. (13). If 
imax = 1, then Umax ~ p 2 /2m and the integral can be taken analytically: 

25 p 2 25 

HUmax = P 2 /2M) = —^- = —U max . (43) 

6 Zm 6 

In this case the energy taken by GWs is of the order of the kinetic energy of PBH and correspondingly 

Q GW ~ Mn bh v 2 /p B H = v 2 . 

Below we will consider more natural situation when M > mpi t> _1 . Integral (42) in the limit of small 

This expression is accurate within 30% up to ^ max = 1. So in what follows we will use this result as 
I(ou max ) ~ 25u max /3, keeping in mind that normally u max = l/r g <C p 2 /2M. 

The fraction of the cosmological energy density of the emitted gravitational waves which has been 
produced during time interval t and t + dt, which is smaller than or comparable to the cosmological time 
t± < t < t ev ~ r BH , can be obtained by the integration of equation (41) over ui from to u max taking 
into account that the energy density of GWs goes with the redshift as (1 + z)~ 4 , and the integration 
over cosmological time, t, which is connected with the redshift by the relation 3 

dt = m , 45 

H* (1 + z) [Q B h*(1 + zf + fi„(l + z)f /2 

where H*, VLbh*, and VL„ are respectively the Hubble parameter, the matter density parameter, and the 
radiation density parameter evaluated at cosmological time t* = tbh, just before the PBH decay. Recall 
that we use the instant decay approximation, so the Universe at t = tbh was still at MD stage. In this 
case all quantities such as and p c are taken at this stage: = 2/3t*, p c = m 2 pl /QnT BH) VL B h* = 1) 
and f2„ = 0. 

We need to calculate the energy density of GWs at the moment of the PBH evaporation. The rate 
of GW production is given by eq. (41). To take into account the redshift of the energy density of the 
gravitational waves we have to divide dpcw/du by (1 + z) 4 , to substitute u — (1 + z)u*, where cu* is 
the GW frequency at t = t B h-, and to express time through the redshift as dt = (3/2)tbh(1 + z)~ 5 ^ 2 dz. 
As a result we obtain at t* = t B h'- 

d PGW (TB H ) = 3 ™] Vrel [p { BH ter) ?TBH(l + z)^ 2 f[ Um (z + l)]d[(l + z)u,]dz ■ (46) 

bm pl 

Here p p l ^ ster ^ is the energy density of the PBHs in the cluster (which is denoted above as pb). Note that 
^(duster) _ cong £ De f ore the PBH decay. We parametrize this quantity as p^ s<er ' ) = Pb^ h (t B h)A(tbh) , 
where P^b H {tbh) = m pJ {^ t bh) ls the average cosmological energy density of PBH and A(tbh) is given 
by eq. (20), see also the discussion above this equation. Function f(oo) is the function of £ = 2mcu/p 2 
in the square brackets of eq. (40). 

To find the cosmological energy fraction of GWs at t = t B h we need to integrate the expression 
above over frequency, using eq. (43), and over redshift and to divide it by the total average cosmological 



3 In this paper we consider flat space with curvature k = and neglect cosmological constant, A = 0. 
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energy density Pbh( t bh) — m pJ {^>' kt bh)- Since we have to average over the whole cosmological volume, 
one factor A disappears and we remain with the first power of A. So the cosmological energy fraction 
of GWs would be: 

V rei \ f A \ f N eff \ (U Tl 



M*- r BH ) » 160 ^fj {w){w) Kit) ■ < 47 > 

Here coefficient Q reflects the uncertainty in the cross-section due to the unaccounted for Sommerfeld 
enhancement [2E1EE]. Note that A may be considerably larger than 10 5 . 

With v re i = 0.1, A = 10 5 , Q = 100, and f max = r" 1 the fraction of the cosmological energy density 
of the GWs emitted by the bremsstrahlung of gravitons from the PBHs collisions, when the Universe 
age was equal to the life-time of the PBH, could reach: 

n GW (T BH )~3.8-10- 17 . (48) 

It looks that for very light PBH, M < 50mpi , the fraction of GW might exceed unity, which is 
evidently a senseless result. However, one should remember the lower bound on the PBH mass (16) and 
that mpi/M < fi p /20 and m P[ jM < 10™ 7 (fi p / 10 ~ 6 )- 

It may be interesting to calculate the contribution to £Igw(t~bh) from the earlier period before the 
cluster formation. The mass density of PBHs at that stage was equal to the cosmological energy density 
but since it was quite high and the effect is proportional to the density squared, the contribution from 
this period might be non- negligible. The result can be obtained from eq. (46), where pbh is taken equal 
to the average cosmological energy density. Since p c evolves with time we need to insert into the integral 
over dz the factor (1 + z) 6 where the redshift is taken from some initial time, presumably ti = t eq , down 
to the moment of the cluster formation, t±. So the energy density of gravitational waves produced by 
bremsstrahlung from t = t eq (9) till t = t\ (15) would be: 

= 32 ^ rel [PBUti))%(l + z)- 1/2 fMz+ l)]d[(l +z)«,,]dz, (49) 

where pg H = m 2 pl /(67rtl) and (1 + z) runs from 1 up to (ti/t eq ) 2 ^ 3 . We have introduced an upper index 
(1) to indicate that this is the energy density of GWs generated before the cluster formation time t = t\. 
The integration over z gives the enhancement factor (1 + z max ) l l 2 = (ti/te^) 1 / 3 . According to eqs. (9) 
and (15), this ratio is A" 1 / 2 ~ 10 2 . Another enhancement factor comes from a larger cosmological energy 
density p^(ti) = p^ (t bh )(t bh /£i) 2 . The other factor /?£#(£i) disappears in the ratio Vt GW = pcw/p^- 
On the other hand, VLqw is redshifted by (tbh Ai) 2 ^ 3 Correspondingly 

n(i) i \ 1 1 a -i/ 2 C 1 ) / \ 1/3 

^gw( t bh) = HA iw / v±l ( tbh\ ^ 
Qgw(tbh) A(t B h) v rel \ t x J 

where the coefficient 11 came from the ratio of the integrals over z of eqs. (46) and (49) and 

The ratio of relative velocities of PBHs before and after the cluster formation, v^/v re i, is tiny, according 
to the estimates of sec. 3, and this introduces another strong suppression factor to the production of 
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GWs at an earlier stage. In accordance with eq. (20) the density contrast rises as A = A(t 1 )(rs_ff/ii) 2 , 
where A(^) is supposed to be large, say, 10 4 — 10 5 due to the fast rise of density perturbations at MD 
stage after they reached unity. Thus the generation of GWs in high density PBH clusters is much more 
efficient than at the earlier stage. 

The density parameter of the gravitational waves at the present time is related to cosmological time 
t% as: 

JW(f ) = (^|4) (fj) , (52) 

where H = 100/i km/s/Mpc is the Hubble parameter and h = 0.74 ± 0.04 [SUl EI]. 
Using expression for redshift (26) and taking the emission time = tbh we obtain: 

/ inn \ 1/3 

Q GW (t ) = 1.67 x 10-% 2 — — n GW (T BH ) . (53) 

\gs{T{TBH))J 

Now using both equations (48) and (53) we find that the total density parameter of gravitational waves 
integrated up to the maximum frequency is: 

where K is the numerical coefficient: 

(£) 

Presumably K is of order unity but since A may be much larger than one, see eq. (20), K may also be 
large. 



0.C10- 21 ir(^) ■ (r,J) 



NejA m ( loo y^ 



100 ; Vl00; \g s {T{r BH ))J 



5 GW from PBH scattering. Classical treatment. 

Classical radiation of gravitational waves by non-relativistic masses is well described in quadrupole 
approximation, see e.g. books [33j [MJ [35] . However, as we have seen, in high density clusters of PBH, 
their relative velocity could be high, see eq. (23), and relativistic corrections may be non negligible. This 
problem was studied by Peters who considered emission of the GWs by two bodies with masses 
M and m, where the former is supposed to be heavy and at rest and the latter, lighter one, moves 
with velocity v. For non-relativistic motion, when t; < 1, and the minimal distance between the bodies 
is larger than their gravitational radii, the energy of gravitational waves emitted in a single scattering 
process is equal to: 

37vr M 2 m 2 v 

IEa " = B<<1 ' (56) 

where b is the impact parameter. 

For the relativistic motion, 1 — v 2 < 1, the emitted energy is: 

M 2 m 2 

6E ° W = 63^(1-^)3/2 • ( 57 ) 

The frequency of the emitted gravitational waves in this process is peaked near u ~ 2ir/St, where St 
is the transition time which, for non-relativistic motion is St = b/v according to ref. [36], while for the 
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relativistic one it is equal to 8t ~ 6(1 — v 2 ) 1 / 2 . For an order of magnitude estimate let us take M ~ m, 
then the radiated energy, as a function of frequency, would be: 

M 4 

^ GIf (w)^— w 3 . (58) 
m pl 

This and the previous equations are true for sufficiently large impact parameter, b ^> r g for which the 
space-time between the scattered PBHs may be considered as flat and their gravitational mass defect 
can be neglected. The energy loss in a single scattering event cannot be larger than 

5E max = , (59) 

where p = Mv re \ is the relative momentum of two scattered PBHs and q is the momentum transfer which 
by an order of magnitude is q — 1/6. Here and in what follows we use non- relativistic approximation. 
So equations (56) and (57) can be true only for 



, , 37n M 2 . . 

b>bmm = yi5^- (60) 

For smaller impact parameters the radiation of gravitational waves would be considerably stronger but 
the approximation used becomes invalid. For the (near) "head-on" collision of black holes a bound state 
of two BH (a binary) or a larger black hole could be formed and the energy loss might be comparable 
to the BH mass due to gravitational mass defect. However, we are interested in gravitational waves at 
the low frequency part of the spectrum, such that they could be registered by existing or not-so-distant- 
future GW detectors. For such low frequency gravitational waves the approximation used here is an 
adequate one. 

The differential cross-section of the gravitational scattering of two PBHs in non-relativistic regime, 
q 2 2M 2 , can be taken as: 

da = — k %- = — 5— bdb. (61) 

q 4 m z pi 

The differential energy density of GWs emitted at time and frequency intervals [t, t + dt] and [u, uj + du] 
respectively can be calculated as follows. The rate of the energy emission by GWs is 

dpcw = dan 2 BH v rel SE GW , (62) 

where we take for SE non-relativistic expression (56). We assume that the impact parameter is related 
to the radiated frequency as u> = 2nv re i/b, as is discussed below eq. (57). So bdb = b 3 dcu/(2iiv re i). So 
we find: 

7AuV re l 2 M 4 du 

d Pow = ^^PBH^r i ^dt. (63) 

The energy density parameter of GW at the moment of BH evaporation can be obtained integrating 
this expression over time and frequency. Thus we obtain: 

«-^0 2 (§)(^)(^). 
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Figure 1: Log-log plot of density parameter today, h^Qcw, as a function of expected frequency today in 
classical approximation for N e ff ~ 100, g$ ~ 100, A& ~ 10 5 , and v re i ~ 0.1 for different values of PBH 
mass M ~ 1 g (solid line) and M ~ 10 5 g (dashed line). 



If we do not confine ourselves to the impact parameter bounded by condition (60) and allow for b ~ r g , 
the energy density of GWs at the moment of PBHs evaporation might be comparable to unity. 

Let us now take into account the redshift of GWs emitted at different moments during the the life- 
time of the high density clusters. The energy density of GWs emitted at some time t is redshifted to 
the moment of BH decay as l/(z + l) 4 . The frequency of GW is redshifted as u = (z + l)w*, where is 
the frequency of GWs at t = tbh- Integration over time or redshift is trivial and we find from equation 
(63) that the energy density parameter of gravitational waves per logarithmic interval of frequency or 
the spectral density parameter, which is defined according to ref. [32] as: 



fW(/;^-^, (65) 

Pc din/ 



at time t = tbh is equal to: 



G ,,/,; TBH ),8,(^)(^)(^)(4)/,. (66) 

Now using equations (31) and (53) we can calculate the relative energy density of GWs per logarithmic 
frequency at the present time: 



h&Gwif; to) « 1-23 • 10- 12 a' 



1 n>i \ 1/2 

10 5 g x 



(67) 
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where a' is the coefficient at least of order of unity: 



/tW\ ( M (NefiV' 2 ( 100_ • 1 : 



«'=lin){w){-m) ) ■ (68) 

It may be much larger if A& 3> 10 5 . 

As we mentioned above, the classical approximation is valid if the impact parameter is bounded from 
below by equation (60). Since the frequency of the radiated GWs is of the order of v/b, the maximum 
present day frequency of GWs, produced at cosmological time t = r BH , for which the classical non- 
relativistic approximation is still valid, would be: 

~ 9 • 10 Hz (— ) ^-^—^j . (69) 

For M = 10 5 g the minimum impact parameter is b m i n ~ 10 13 cm. The frequency of the order of 1 
Hz today corresponds to the impact parameter 6 orders of magnitude larger. If we demand that the 
impact parameter should be smaller than the average distance between PBHs in the clusters, then using 
equations (22) and (60) we find that it can be true if the following condition is fulfilled: 

6 Energy loss of PBHs 

We calculate here the total energy loss of PBHs in the high density clusters, in order to understand 
how probable could be the formation of the PBH binaries. First, let us estimate the total energy loss 
of PBHs due to the graviton bremsstrahlung. The loss of the kinetic energy per unit time due to the 
graviton emission is: 



/ dE kin \ r ( da\ 

—7— = n B HVrel I dbJ U [ — ] 

V at J brem n \ aUJ / brem 



'max / -i \ 



where io max is defined in sec. 4, The total loss of kinetic energy of a single PBH during the time interval 
equal to the PBH life-time, 5E)~ in = —E kin TBHi normalized to the original kinetic energy of the PBH 
can be estimated as 

SEkin „ . „ 4 frnpi\ 2 



Ekin 

where 



6.10%(^) 2 , (72) 



Clearly the energy loss is essential for very light PBHs which could form dense clusters only if Q p is 
sufficiently high, see eq. (16). 

The energy loss due to classical GW emission might be somewhat more efficient. According to the 
previous section the energy loss by a single PBH per unit time is: 

AE class = n BH v f db[^-) 8E(b) , (74) 

i. \ a0 / class 
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where 5E(b) and b min are given respectively by eqs. (56) and (60). 

Taking the integral over b and time we find for the fractional energy loss of PBH due to classical 
emission of the gravitational waves: 



AE, 



class 



E, 



0.9-10° 



kin 



A 6 Neff mpi 
10 5 100 M 



(75) 



One should remember however that this energy loss comes from the PBHs scattering with rather large 
impact parameter b > b min . For smaller b, when the simple approximation used in this work is inappli- 
cable, the energy loss might be much larger. Moreover, according to eqs. (9), (15), and (20) the density 
amplification factor A& may be much larger than 10 5 : 



A 6 (t B h) = 10 4 A(^ 



A ft 4 

in p 



100 



N, 



eff 



M 

mpi 



(76) 



where we may expect e.g. that A(£i) ~ 10 5 , A in ~ 10~ 4 , and Q p ~ 10~ 6 . 

PBHs in the high density clouds could also loose their energy by dynamical friction, see e.g. book 
A particle moving in the cloud of other particles would transfer its energy to these particles due to their 
gravitational interaction. However, one should keep in mind that the case of dynamical friction is 
essentially different from the energy loss due to gravitational radiation. In the latter case the energy 
leaks out of the system cooling it down, while dynamical friction does not change the total energy of the 
cluster. Nevertheless a particular pair of black holes moving toward each other with acceleration may 
transmit their energy to the rest of the system and became gravitationally captured forming a binary. 

For an order of magnitude estimate we will use the Chandrasekhar's formula which is valid for a 
heavy particle moving in the gas of lighter particles having the Maxwellian velocity distribution with 
dispersion a. The deceleration of a BH moving at velocity vbh with respect to the rest frame of the gas 
is given by 



d _ 
dt VBH 



—An G N Mbh Pb In A 



vbh 



'BH 



erf(X) 



2Xexp(-X^ 



where X = Vbh/(v2ct), erf is the error function, p b is the density of the background particles, 
In A m ln(M*/MB#) is the Coulomb logarithm, which is defined as 



(77) 



and 



In A = In 



c m Pl a 



M B h + m 



Here b max is the maximum impact parameter, a 2 is the mean square velocity of the gas and m is the 
mass of particles in the gas. Numerical simulations show that b max can be assumed to be of the order 
the radius of the cloud, Rb, which is given by equation (21). Since a 2 ~ Mb / \m 2 Pl Rb) , a reasonable 
estimate of A is Mb/ Mbh- 

Equation (77) was solved in ref. [38] in two limits v > a and v < a. In both cases the characteristic 
dynamical friction time was of the order of: 



T~DF 



S 4 

o m Pl 



An M B h pb In A 



a 

01 



25 



_ln(10- 6 /fi p ) 



100 



N. 



eff 



M\ /10 6> 



(7f 



For PBH masses below a few grams dynamical friction would be an efficient mechanism of PBH cooling 
leading to frequent binary formation. Moreover, dynamical friction could result in the collapse of small 
PBHs into much larger BH with the mass of the order of M& (18). This process would be accompanied 
by a burst of GW emission. 
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7 Gravitational waves from PBH binaries 



Binary systems of PBH could be formed with non-negligible probability in the high density clusters. As 
we have seen in the previous section, PBHs could loose their energy due to emission of gravitational 
waves and due to dynamical friction [37] . As a result they would be mutually captured. Determination 
of the capture probability is a complicated task, which could probably be solved by numerical simulation. 
Since it is outside of the scope of the present work, we simply assume that the mass or number fraction 
of PBH binaries in the high density bunches of PBH is equal to e, where e is a dimensionless parameter 
which is hopefully not too small in comparison with unity. 

Gravitationally bound systems of two massive bodies in circular orbit are known to emit gravitational 
waves with stationary rate and fixed frequency which is twice the rotation frequency of the orbit. In 
this approximation orbital frequency, w rf>, and orbit radius, R, are fixed. Luminosity of GW radiation 
from a single binary in the stationary approximation is well known, see e.g. book |33j : 

. 32M 1 2 Mf(M 1 +M 2 ) 32 2 fM c u orb \ w/3 

where M\, M 2 are the masses of two bodies in the binary system and M c is the chirp mass which is 
defined as 

(Mi M 2 ) 3/5 , v 



and 



orb 



M l + M 2 
m 2 Pl R 3 



In the case of elliptic orbit with large semi-axis a and eccentricity e the luminosity is somewhat larger 
(if R = a): 

32M 1 2 M 2 2 (M 1 + M 2 ) t 73e 2 + 37e 4 \ 



5a 5 m 8 Pl (1 - e 2 ) 7/2 V 24 96 



The emission of GWs costs energy which is provided by the sum of the kinetic and potential energy of 
the system. To compensate the energy loss the radius of the binary system decreases and the frequency 
rises making the stationary approximation invalid. As a result the system goes into the so called inspiral 
regime. Ultimately the two rotating bodies coalesce and produce a burst of gravitational waves. To 
reach this stage the characteristic time of the coalescence should be shorter than the life-time of the 
system. In our case it is the life-time of PBH with respect to the evaporation. 

In the inspiral regime the initially circular orbit may remain approximately circular if radial velocity 
of the orbit, R, is much smaller than the tangential velocity, u or bR. This regime is called quasi-circular 
motion and is valid as long as (see e.g. book [39]): 

^orb < ^Irb- ( 83 ) 

Equation (83) can be translated into the lower bound on the radius of the orbit: 



m 2 Pl 



*4) 
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which is the condition of the validity of the Newtonian approximation. It was shown by Peters [40J that 
the orbits with initial e$ = would remain quasi-circular as far as condition (83) is fulfilled, while for 
the orbits with eo 7^ the eccentricity rapidly approaches zero due to back reaction of the gravitational 
radiation. 

Most probably binaries are formed in elliptic orbits with high eccentricity. However in the calculation 
of the GW emission by binaries we assume for simplicity that all orbits are circular. The result would 
be a lower bound on GW emission, hopefully not too far from the real case. 

In what follows we will consider both stationary and inspiral regimes since they both might be 
realized for different values of the parameters. We will use the instant decay approximation, when the 
PBH mass is supposed to be constant till t = tbh and then BH would instantly disappear. The case of 
the realistic decrease of PBH mass will be considered elsewhere. 

The stationary orbit approximation would be valid if time of coalescence, r co , would be much larger 
than the BH life-time, r co > tbh- The former can be found as follows (see e.g. book [33]). According 
to the virial theorem the total (kinetic plus potential) energy of the system is £ = — MiM 2 /(2i?m 
Since luminosity (79) is L s = —dS/dt, the radius varies with time according to 



2 \ 
pu 



Correspondingly 



. 64M 1 M 2 (M 1 + M 2 ) 

R= 5i?%4 • ( 85 ) 



/?</) n it [ to + Tco * V /4 , (86) 



where Rq is the initial value of the radius, t is the initial time, and the coalescence time is given by: 

(87) 



256MiM 2 (Mi + M 2 ) 

The condition r co > tbh can be translated into the lower bound on R (for Mi = M 2 ): 

fl>iU. = «.10»g) I "(^ 5 ) V, r.. (88) 

Keeping in mind that the frequency of GWs emitted at circular motion of the binary is twice the orbital 
frequency, f s = u or b/n we find from equation (81) that lower bound (88) leads to the following upper 
bound on the GW frequency: 

/. < cWtt « 2 ■ 10 24 Hz ( ^ ( ^) V4 . (89) 



^ 100 y \ m 

On the other hand, the radius of the binary orbit should be smaller than the average distance between 
PBHs in the cluster (22) and probably quite close to it. Using eqs. (22) and (88) we find: 

So it seems natural that R m in "C db and the PBH binaries should be mostly in the quasi- stationary 
regime. R m in would be equal to db roughly speaking for quite large mass fraction of the produced PBHs, 

a p > 10- 3 . 
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The condition R r , 



db gives a lower bound on orbital frequency, oo or b- 

1/2 / A x 3/2 / ^ \ 2 



CU or 6 > CO 
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9.4- 10 17 sec -1 
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io- 



10 5 g 
M 



(91) 



During the inspiral phase, for which r co < tbh, we expect that binaries emit GWs in the frequency 
range: 



2 • 10 24 Hz 



N, 



eff 



100 



3/8 



10 5 g 
M 



7/4 



< / < 0.6 • 10 33 Hz 



10 5 g 
M 



(92) 



The upper bound corresponds to u; ~ l/r g . 

The frequency spectrum of the gravitational waves in inspiral but quasi-circular motion can be found 
in the adiabatic approximation as follows. Since the gravitational waves are emitted in a narrow band 
near twice the orbital frequency, the spectrum of the luminosity (79) can be approximated as: 



dE 



32M 1 2 Mf(M 1 + M 2 ) 



5R 5 (t)m 



S (uj — 2u or b(R)) dco 



(93) 



pi 



To find the energy spectrum we have to integrate this expression over time from initial time, t m i n = to, 
to maximum time t max = min[rBH + t p , r co + to], where to an d t p are respectively the time of the binary 
formation (it may be different for different binaries but here we neglect this possible spread) and the 
time of PBH formation (it is different for PBH with different masses). Note that the coalescence time, 
t co is also different for binaries with different initial radius Rq. 

Using eqs. (81) and (85) and the expression dt = (d-R/dt) _1 (d-R/do; or b)da; or fe, we find: 



dE 
din a; 



2 V3 w 2/3 



M X M 2 



m^ 3 (M 1 + M 2 )V3 



(94) 



in agreement with refs. [39], |4T]. This expression is valid for the frequencies in the interval determined 
by eq. (81) with R max = R and R min = R{t max ). 

In expression (94) we have not taken in account the redshift which is different for different frequencies 
and thus this leads to spectrum distortion. According to eqs. (81) and (86) frequency u is emitted at 
the time moment: 



t(cu) = t + 



where 



1 



M 1 + M 2 



OJr. 
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pi 
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8/3' 



(95) 



R 



-3/2 



(96) 



is the minimal frequency emitted at initial moment t = t . To the moment of the PBH evaporation the 
frequency of the GWs emitted at t = t(u) is redshifted by the frequency dependent factor: 



00 



CO* 



l + z(uo) 



1 2/3 



t p + tbh 



lo, 



(97) 



where co* is the frequency of GWs at t — t p + r BH . This equation implicitly determines oo as a function 
of oo*. 
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The spectrum of the gravitational waves at PBH evaporation can be obtained from eq. (94) dividing 
it by (l+z) (the redshift of the graviton energy, E) and with substitution u> = (z+1)lu*. Correspondingly 

1 — uj*{dz/du) 

As a result we find: 

dE* 2 1 /3 W 2/ 3 Ml M 2 [l-u^dz/du)]' 1 



dlnw, 3 m J/ 3 (M 1 + M 2 )V3 (l + z)V3 



(99) 



Here z(oj) should be taken as a function of according to eq. (97) and varies between u min and 
^max divided by the corresponding red-shift factor. In particular, uj*( m i n ) = to m i n [to/(t p + r BH )} 2 ^. Note 
that Rq enters explicitly into eq. (99), while in eq. (94) it enters only through the limits in which oo 
varies. Because of that the frequency spectrum depends upon the distribution of binaries over their 
initial radius, Rq. As is shown below, it is especially profound in the case of long coalescence time when 
the frequency spectrum of a single binary with fixed R is close to delta-function. 

In the stationary approximation, when the change of the orbit radius can be neglected, we expect that 
a single binary emits G Ws in a narrow band of frequencies close to twice the orbital frequency. However 
the distribution of binaries over their initial radius, driBiN = F(Ro)dRo spreads up the spectrum. Here 
driBiN is the number density of binaries with the radius in the interval [R ,R + dR ]. Since in this 
approximation the radius is approximately constant we do not distinguish between R and R . The 
cosmological energy density of the gravitational waves emitted per unit time is equal to: 

^ = 2^=™*^., (100) 

where n b BH is the number density of PBH in the high density bunch (cluster), n c BH is the average 
cosmological number density of PBH, R = R(uj or b) according to eq. (81), and we used the relation 
dR = —2(R/3) (dui/ui). Distribution, F(R), is normalized as: 



/ 



dRF(R) = n BIN = en b BH . (101) 



We assume for simplicity that F(R) does not depend upon R in some interval [Ri, Rq\ and vanishes 
outside it. So F(R) = en b BH /(R 1 - R 2 ). 

A more realistic fit to the PBH distribution over radius could be a Gaussian one: 

F{R) = —^en BH exp [-(R - (R)) 2 /2a 2 ] , (102) 
v 2tt a 



where a is the mean-square deviation of R from the average value (R). 

The small factor n c BH /n BH enters eq. (100) because we are interested in the cosmological energy 
density of GWs averaged over the whole universe volume. The cosmological number density of PBH is 
expressed through their energy density as n BB = p B n/M = p c (t)/M. The number density of binaries 
in the cluster is parametrized according to: 

n mN (t) = e(t) n b BH (t) = e(t) p c (t)A(t)/M , (103) 

where, we remind, p c (t) is the total cosmological energy density and A(i) = pb/p c ^> 1 is the density 
contrast of the cluster. The time dependence of n b BH disappears when the cluster reaches the stationary 
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state, see discussion in Sec. 2, and A(t) evolves according to equation (20). When the stationary orbit 
approximation is valid, e remains constant. 

Collecting all the factors and integrating eq. (100) over time with an account of the frequency 
redshift, u = u;*(l + z) and the total redshift of the energy density of GWs, paw{t*) = PGw(t)/(l + z) 4 , 
we find: 



i (stat) 
d PGW 1 



2 7/3 



n BH {r BH ) 



n 



BH 



(M?Mp(T BH + t p ) 



F{R)u 5 J 3 du* J rr n / 6 dx, 



(104) 



where x = a(t)/a(t*) = 1/(1 + z), x mm = a(to)/a(t*), to is the time moment of binary formation and we 
make use of equation (45). Dividing this result by the critical energy density just before PBHs complete 
evaporation, ubh{jbh) ~ Pc( t bh)/M, we find the cosmological fraction of the energy density of GWs 
at t — t B h per logarithmic interval of frequency / = u/(2n) (below we assume that all BHs have equal 
masses, M): 



3 ■ 2 17 / 3 e-(t„ + t bh ) (irf*M\ 8/3 
R,-R 2 \m 2 Pl ) [ 
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(105) 



where for the sake of a simple estimate we assumed that F(R) = const. We assume also that all the 
binaries are formed at the same time, t <C t B h and so x min <^ 1. Note that the frequency of GWs 
coming from the binaries with radii between R\ and R2 is confined according to eq. (81). 

To make an order of magnitude estimate of the fraction of the energy density of GWs at the moment 
of PBH evaporation we take (i?i — R 2 ) ~ R\ ~ R(w), where R(cu) is determined by equation (81) and 
take into account that the stationary approximation is valid if the radii of the binaries are bounded 
from below by eq. (88). Hence, if the stationary regime is realized, the spectral density parameter today 
would be: 
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(106) 



The expected range of the present day frequencies of the GWs from the binaries in the stationary 
approximation is given by eqs. (91) and (89). The emitted frequency is determined by the binary radius, 
so a single binary emits GWs with a very narrow spectrum. However, the distribution of binaries over 
their radius could lead to a significant spread of the spectrum. In principle the frequencies emitted may 
have any value in the specified above range. The minimal present day frequency of such GWs today 
can be found by plugging eq. (91) into eq. (31): 
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(108) 



For binaries formed with R > R m i n , see equations (31), (88) and (89), the frequency of emitted GWs 
today is bounded from above by: 
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Figure 2: Log- log plot of density parameter today ti^VLcw as a function of expected frequency today for 
PBHs binaries in the stationary approximation for (3 ~ 1, e ~ 10~ 5 , N e ff ~ 100, g$ ~ 100, PBH mass 
M ~ 10 7 g (solid line) and M ~ 1 g (dashed line). 



Let us estimate now the energy density of GWs in the inspiral case, when r co < tbh and the GW 
emission from a single binary proceeds in a wide range of frequencies due to shrinking of the binary 
radius. The radiation frequency spans from f Si min, which is the GW frequency at the initial PBH 
separation, to f s ,max which corresponds to GWs emitted at R ~ r g . The energy spectrum of GWs is 
given by eq. (94) where, in what follows, we change to cyclic frequency, / = u/2tt. 

After the cluster evolution was over, the number density of PBHs in high density clusters remained 
approximately constant till the PBH evaporation, but in the inspiral phase the fraction of binaries, e(t), 
decreased due to their coalescence. So the tail of the distribution function at small initial Rq is eaten up, 
and the average value of R drops down. In distribution function, F(R ), we have to substitute instead 
Rq its expression through R and time according to 



Rq -> 



R' + 



^256M 1 M 2 (M 1 + M 2 



5mp ; 



(t - 1 ) 



■,1/4 



(110) 



with the corresponding change of R^dRo — > R dR. 

To calculate the cosmological energy fraction of GWs at the PBH evaporation moment we can proceed 
along the same lines as we have done deriving eq. (99) introducing additional factor F(R )dRQ which 
depends upon time according to eq. (110). However, at the level of calculations in the present model 
with many unknown parameters it can be sufficient to neglect such subtleties and to use a simplified 
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estimate: 

where e co is the fraction of binaries with coalescence time shorter or equal to PBH life-time. For an 
estimate by an order of magnitude we assume also that the number of binaries is independent on the 
redshift. To some extend the decrease of the binary number may be compensated by their continuous 
formation. We neglect possible difference of binary masses and take Mi = M 2 . We approximately take 
the redshift into account from the moment of the coalescence to the PBH decay, (z co + l) « (tbh/ t co) 2 ^ 3 - 
This corresponds to the assumption that the binaries radiated all GWs only at the moment of r co . So 
the /* = /(l + z co ). Thus we obtain as an order of magnitude estimate: 

fW(/., t bh ) = "-f (^-) (z co + I)" 1 /* . (112) 
3 V m pi J 

Using equations (31) and (53), we find that the energy density parameter of gravitational waves 
today is equal to: 

where we neglected possibly weak redshift dilution of GWs by the factor (t co /t bh ) 2 / 9 . 

If the system goes to the inspiral phase, then according to equation (92) we would expect today a 
continuous spectrum in the range from f min ~ 0.9 • 10 7 Hz to f max ~ 3 • 10 14 Hz. However if we take into 
account the redshift of the early formed binaries from the moment of their formation to the PBH decay, 
the lower value of the frequency may move to about 1 Hz. 



8 Gravitons from PBH evaporation 

In the previous sections we have considered only gravitational waves emitted through mutual acceleration 
of PBHs in the high density clusters. On the other hand PBHs could directly produce gravitons by 
evaporation. This process in connection with creation of cosmological background of relic GWs was 
considered in ref. [6] and later in ref. [12]. In the last reference a possible clumping of PBHs at the 
matter dominated stage was also considered. Though such clumping does not influence the probability 
of the GW emission by PBHs, it may change the mass spectrum of PBHs due to their merging. 
The PBHs reduce their mass according to the equation: 

( t-t \ 1/3 

M(t) = M (l-—±) , (114) 
V r BH J 

where Mo is the initial mass of an evaporating BH and t p is the time of BH production after Big Bang. 
Equation (114) shows that the BH mass can be approximately considered as constant till the moment 
of the evaporation and may be approximated as 6(t — t p — tbh)- Due to evaporation a BH emits all 
kind of particles with masses m < Tbh and, in particular, gravitons. The total energy emitted by BH 
per unit time and frequency u (energy) of the emitted particles, is approximately given by the equation 
(see, e. g. book [43J): 

dE \ 2N eff M 2 u, 3 (n5) 



(dtdw J 7r m Pl c^/Tbh - 1 
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where T is the BH temperature (2). Due to the impact of the gravitational field of BH on the propagation 
of the evaporated particles, their spectrum is distorted [S] by the so called grey factor g(oj), but we 
disregard it in what follows. 

Let us now estimate the amount of the gravitational radiation from the graviton evaporation. After 
their production PBHs started to emit thermal gravitons independently on the PBH clustering. Hence 
the thermal graviton emission depends only on PBH number density, n BB - The energy density of 
gravitons in logarithmic frequency band emitted in the time interval t and t + dt is 

dp GW (u;t) 2 ( dE \ 

-10 n BH {t) \-r— )dt, (116) 



do; \dtduj 

where factor 10~ 2 takes into account that about one percent of the emitted energy goes into gravitons. 
The density parameter of GWs per logarithmic frequency interval at cosmological time = tbh can be 
obtained by integrating expression (116) over redshift with an account of the drop-off of the graviton 
energy density by (1 + z)~ 4 and the redshift of the emitted frequency so that at £* = r BH : uj = u*(l + z). 
Note that in the instant decay approximation the BH temperature remains constant. One has also to 
take into account that the number density of PBH behaves as riBH{t) = n p (t p )(l + z) 3 , so if we normalize 
our result to fi BB \{r bh)i the integrand should be multiplied by (1 + z) 3 . Finally we obtain: 

J] = 4 ( 3T BH) Pbh{t BH ) I 7fT~ , 117 

din to* nm^t \T BH J 



where 



and 



uj* \ = Zm r d^(i + ^) 1/2 

Tbh) ~~ i exp [{z + 1)ujJT bh } - 1 ' 



us) 



where the effective time of integration is equal to 3t bh because of the instant decay approximation. One 
can check that in this case the total evaporated energy would be equal to the PBH mass. 
The spectral density parameter of GWs at t — t B h is equal to: 



2.9 ■ 10 3 MV 



uj* 



n GW (uj*; r BH ) « 3 *- 1 I — — 1 . (120) 

ixm pi \T BH J 

The spectrum is not a thermal one, though rather similar to it. It has more power at small frequencies 
due to redshift of higher frequencies into lower band and less power at high u*. The spectral density 
parameter reaches maximum at ujP eak /T BH = 2.8. Accordingly the maximum value of the spectral 
density parameter when PBHs completely evaporated is equal to: 

n p G t(ujr k ;r BH )^3.8-io- 3 . (121) 

Integrating equation (120) first over uj* and then over redshift, we find that the total fraction of 
energy of GWs is 0.006 which is reasonably (in view of the used approximations) close to the expected 
0.01. At BBN the energy fraction of such GWs would be about 0.005. So the total number of additional 
effective neutrino species would be close to 0.045, where 0.03 comes from neutrino heating by e + e~ 
annihilation and 0.01 comes from the plasma corrections (see e.g. review [44J). Of course the GWs 
produced by the considered mechanism are safely below the BBN bound [45]. Using equation (53) 
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Figure 3: Log-log plot of the density parameter per logarithmic frequency, hlQ GW (f; t ), as a function 
of frequency today, /, for the case gs ~ 100, N e fj ~ 100, black hole mass M = 1 g (solid line) and black 
hole mass M = 10 5 g (dashed line). We can see that the spectrum has a maximum which is sharp and 
of order hlVt GW (f peak ) ~ 10~ 7 . 

and taking into account the redshift from t = tbh to the present time, we find that the total density 
parameter of GWs today due to PBH evaporation would be about 10~ 7 . 

The total energy density of GWs from the PBH evaporation is quite large but it is concentrated at 
high frequencies. According to eq. (31) the redshifted peak frequency emitted at time = tbh becomes 
today: 

= 2 . 10 » Hz fecw) 1/12 (m\ 1/2 («) 1/2 . (122) 

V 100 J \N e „J \WgJ ^ 1 

The energy density of GWs at small / drops down in accordance with equation (121). The spectral 
density today can be calculated from equation (120) with an account of the redshift to the present day: 

«^,3 6 ,0^(^) 2 (^) 2 (^) 4 ,(^), m 
where we used u = 2irf and T is the BH temperature redshifted to the present time: 

To = Tbh = 4 ,3 . ^ H Z («» ) 1/12 ( ™ ) (") " 2 . (1 24) 

. a(*o) J \9s{T{tbh))J \NeffJ \10 5 gJ 
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9 Stochastic background of gravitational waves. An overview 



Stochastic background of relic gravitational waves can be produced by several mechanisms. The theo- 
retical predictions are model depended due to the uncertainties in the cosmological framework and on 
the values of the redshift from the production epoch. Below we briefly describe some of the production 
scenarios. For a more detailed review on stochastic background of GWs production mechanisms and 
their spectra the reader can consult more specific ref. [SI H7J Pf8] . 

Inflationary models. It was established long ago that gravitational waves could be produced in 
cosmology due to an amplification of vacuum fluctuations by external gravitational field (quantum 
particle production). It was first studied by Grishchuck [3 J and first applied to an inflationary model by 
Starobinsky[3]. The gravitational waves could be quite efficiently produced at inflation. Their spectrum 
at large wavelengths is independent on the details of inflationary models. The frequency band of these 
gravitons today is quite wide and the associated density parameter is very low. The predicted density 
parameter of gravitational waves in the frequency range from 3 x 10~ 18 Hz< / < 10 -16 Hz is 

^.(/(^.n.lO-^) 2 ^) 2 , (125) 

while in the frequency range 2 x 1CT 15 Hz< / < f m ax — 10 9 Hz the spectrum is flat and the density 
parameter is 

h 2 n GW (f) ~ 6.71 • 10- 14 (j^^) , (126) 

where H is the Hubble parameter at inflation. 

A near scale-invariant spectrum over a wide range of frequencies is a key prediction of the standard 
inflationary model (HJ [50] . The relative amplitude of GWs spectrum to density perturbations spectrum 
is usually expressed in terms of the ratio, r, of tensor to scalar perturbations. From observations of 
WMAP, the current limit on B-mode of the CMB polarization demands r < 0.22 which rule out some 
models of inflation [511 [52] . The spectrum of GWs can be expressed in terms of the tensorial spectral 
index, n t , and is almost flat in the frequency range 2 x 10~ 15 Hz < / < f m ax — 10 10 Hz. The density 
parameter is proportional to a power of the frequency: 

h 2 Q GW (f) oc f n \ (127) 

Since the tensorial spectral index is negative, n t < 0, the spectrum is decreasing rather than flat. 
Depending on inflationary model the value of the tensorial spectral index changes and there are some 
models which predict r ~ 10~ 3 . 

Pre-heating phase at the end of inflation. At this stage the energy of scalar field <fi is spent to generate 
new particles and heat the Universe. The first estimate of the density parameter of GWs during the pre- 
heating phase was done by Klebnihkov and Tkachev |53j who found the density parameter of the order 
hy^cw ~ 10~ n for the gravitational waves with the present day frequency / ~ 10 6 Hz, in the models 
with quartic potential, A0 4 . Later, this mechanism was reconsidered by Easther and Lim [5^1 |5"5] who 
studied the models with the potentials of the form A0 4 and m 2 2 . The authors have found numerically 
that h^Qcw ~ 10~ 10 in the frequency range / ~ 10 8 — 10 9 Hz. 

First order phase transitions. At the end of inflation, first-order phase transitions could have 
generated a large amount of gravitational waves. At such transitions the bubble nucleation of true 
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vacuum states and percolation can occur accompanied by the bubble collisions. In a series of pa- 
pers EZl EEl ES] the energy of gravitational waves generated from bubble collisions at strongly 
first-order phase transitions was estimated and the results were later extended to the electroweak first- 
order phase transitions. The amount of GWs from strongly first-order phase transition at its end is of 
the order 1.3 • 10~ 3 (t/H), where r is the duration of the phase transition, H is the Hubble constant, 
and the peak frequency is u^ eak = 3.8/r. The present day density parameter of GWs produced at the 
electroweak first-order phase transition was found to be of the order flow ~ 10~ 22 with characteristic 
frequency / ~ 4 • 1CT 3 . Since later it has been found out, that there is no first order electroweak phase 
transition in the standard model [60J, the mechanism was reconsidered by Grojean and Servant [6T] . 
The authors estimated the GW production in the temperature range 100 GeV-10 7 GeV. The spectrum 
of the GWs today in this temperature range extends from 10~ 3 Hz to 10 2 Hz. The associated density 
parameter was found to be quite large, hlflcwifpeak) ~ 10~ 9 depending on the parameters of the model. 

Topological defects and cosmic strings. Practically in all inflationary models the gravitational wave 
spectrum is almost flat in the frequency range from 10~ 15 Hz< / < f max — 10 10 Hz with some variations 
coming from pre-heating and reheating phases for which the frequency is peaked near GHz region. There 
are other mechanisms of GWs production e.g. by cosmic strings which predict almost flat spectrum in a 
wide range of frequencies. Many of the proposed observational tests for the existence of cosmic strings are 
based on their gravitational interactions [521 E3]- Particularly interesting are GWs produced by closed 
string loops which oscillate in relativistic regime. The spectrum of the gravitational waves produced by 
such relativistic oscillations is almost flat in the region 10 -8 Hz< / < f max — 10 10 Hz with a peak at low 
frequency near / ~ 10 -12 Hz. The density parameter in the frequency range / ^> 10~ 4 Hz, according 
to ref. [61] . is equal to: 

where G/i, a and 7 are respectively the string tension, the initial loop size as a fraction of the Hubble 
radius and the radiation efficiency. From the pulsar timing data the authors of ref. |65j constrained the 
density parameter of GWs from the cosmic strings in the frequency range / > 10~ 6 Hz and put the 
limit 

hln GW {f) < 10- 8 . (129) 

It is generally assumed that at the end of inflation the inflaton oscillates and eventually decays. If 
non-topological solitons, the so called Q-balls, are produced at the inflaton decay, such Q-balls could be 
a source of GWs. According to the calculations of ref. [66] the density parameter of such GWs would 
be of the order of ti^VLcw ~ 10~ 9 with a peak frequency / ~ 10 10 Hz. 



10 Gravitational waves detectors. Present status 

For most of the models mentioned above, the stochastic background of GWs is beyond the sensitivity 
of the current and planned interferometers. We have seen that inflationary models predict almost flat 
spectrum of GWs in a wide range of frequencies. There is a narrow band of frequencies of this background 
that falls into the range of the present detectors such as LIGO and VIRGO. Unfortunately, the density 
parameter predicted by inflationary models is too low to be detected by the present detectors. Almost 
all the models mentioned above predict the density parameter of the order h^flaw % 10~ 5 and actual 
LIGO and VIRGO are not able to detect such a quantity because of the frequency dependence of the 
density parameter. This can be seen from the relation between the expected amplitude of stochastic 
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gravitational waves h c (f) with the density parameter as presented in the ref. [17] 

h c (f) = 1.3 x 10- 18 ^W(7) (jy) • (130) 

Present detectors such as LIGO and VIRGO with enhanced technologies operate in the frequency range 1 
Hz - 10 4 Hz and can reach respectively the strain sensitivity h rms ~ 10~ 23 Hz -1 / 2 and h rms ~ 10~ 22 Hz -1 / 2 
in the frequency band / ~ 10 2 — 10 3 . The planned detectors such as lAdvanced LIGO. [Advanced VTR GQ 




Figure 4: log[/iQ^Gvy(f)] vs - l°g(f P 2 ]) f° r different models of production of stochastic background of 
GWs as given in ref. [45*] . 

and LISA have better chances to detect this stochastic background. In fact, LISA can reach the density 
parameter of the order h^Qcw 1CT 11 at frequency / ~ 1CT 3 Hz and Advanced LIGO can reach a 
h^GW ^ 10~ 9 at frequency / ~ 10 2 Hz. These planned detectors can register the stochastic background 
of GWs coming from cosmic strings and the pre-bing bang stage. The gap between LISA and the ground 
based detectors will be covered by DECIGO/BBO detectors which will operate in the frequency range 
from 0.1 Hz to 10 Hz and have 10 3 better sensitivity than LISA from 0.1 Hz to 1 Hz [67]. DECIGO will be 
able to observe the stochastic background of GWs produced at inflation and can reach IiqVLgw ~ 10~ 20 at 
/ ~ 1 Hz after 3 years of observation [68l EH] . All the above mentioned GWs detectors cover a frequency 
range 10~ 7 Hz - 10 3 Hz and the high-frequency range will hopefully explored by future high-frequency 
GWs detectors. The principle of a high-frequency detector is based on the electromagnetic-gravitational 
resonance first proposed by Braginsky and Mensky [7_3j [7H [75] . Actually there is a renewed interest on 
these new detectors which a prototype has been constructed at Birmingham University [70j (TTJ [72] and 
which reaches a strain sensitivity of the order h rms ~ 10~ 14 Hz^ 1 / 2 at / ~ 10 8 Hz. The main goal of 
this detector is detection of high-frequency stochastic background of GWs from the early Universe and 
black hole interactions in higher dimensional gravitational theories. 
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11 Summary and results 



We have analyzed the formation and evolution of light primordial black holes in the early Universe 
which created a transient matter domination regime in contrast to the present standard cosmology, 
where the early Universe after inflation was normally radiation dominated. PBHs with masses less than 
M ~ 10 8 g evaporated before primordial nucleosynthesis leaving no trace. Thus the fraction of the 
energy density of such PBHs, Q p , in this case is a free parameter of the model, not constrained by any 
existing observations. 

At MD stage the PBHs could form high density clusters which would be efficient sources of the 
primordial GWs. PBHs could have dominated the Universe for a short time of the order of their 
lifetime, t bh , generating relic gravitational waves by various mechanisms of their mutual interactions 
as well as due to their evaporation. In the former case we have shown that production of GWs is most 
efficient after BH density started to dominate over radiation. After that moment, high density clusters 
of PBHs could have been formed, leading to an efficient production of GWs. To survive till cluster 
formation the PBH mass at production must be bounded from below by M ~ 4 • 10~ 5 g Q~ 1 A~^ 4 N^ 

which leads to a lower bound Q p > 10~ 14 A i ^ 4 Nl^. According to the standard cosmology the amplitude 
of primordial density perturbations is of order of ~ 10~ 4 , which in our case leads to a lower bound on 
the density parameter of PBHs, tt p > 10~ n . 

In this context we have calculated the density parameter of GWs today from scattering of PBHs 
in both classical and quantum regime, GWs emission from binaries, and from black hole evaporation. 
We have shown that a substantial amount of gravitational waves has been emitted by all mechanism 
considered here. In the case of scattering of PBHs we considered only scattering between them neglecting 
the possibility of PBHs mergers, which results in an underestimate on ti^fl GW . Even in this case the 
density parameter is substantial at high frequencies reaching values of the order of IiqQgw ~ 10~ 9 at / ~ 
GHz for classical scattering and the total density parameter h^flcw ~ 10~ 10 for very light primordial 
black holes. In the low frequency limit the density parameter in the classical case is of the order of 
h^Vl GW ~ lCT 17 — 1CT 20 in the frequency range / ~ 1CT 1 — 10 2 Hz which falls into the detection band 
of DECIGO/BBO. 

The number of PBHs that form binaries after cluster formation is subjected to uncertainties and in 
this paper we parametrized it through factor e. The exact value of this parameter could be calculated 
elsewhere by numerical calculations. Since the density in such clusters is very high we expect that e is 
not very small in comparison with unity. In fig. 2 the expected value of the density parameter today is 
presented. We can see that a large amount of gravitational waves has been emitted in the high frequency 
regime with IiqQgw ~ 10~ 14 — 10~ 12 at frequency / ~ 10 10 Hz depending on the BH initial mass. In 
the low frequency part of the spectrum the spectral density parameter is utterly negligible making it 
impossible to detect GWs produced by this mechanism at present and probably in the near future. 
In our derivation we have considered both stationary and inspiral phases of binaries leading to a wide 
range of the frequencies emitted. We have considered only binaries in circular orbits and the problem 
with elliptical orbits will be treated later. If elliptical orbits were frequent, the amount of GWs will 
be presumably higher over a wide range of frequencies. We assumed that all binaries are formed with 
initial radius less than the average distance between PBHs and greater than the gravitational radius r g . 
In this case the frequency spectrum has a cutoff in both low and high frequency bands of the spectrum. 

Another mechanism of graviton productions considered here is the PBHs evaporation. This mecha- 
nism is independent on the structure formation during the PBH domination. In fig. 3 we show the density 
parameter as a function of frequency for BH masses 1 g and 10 5 g. Having a near blackbody spectrum, 
the frequency of the emitted gravitons can have any value, but unfortunately the GWs spectrum has 
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a peak in the high frequency region which today make a substantial contribution into the cosmological 
energy density of the order of hlfl GW (f peak ) ~ 1CT 7 . 

The mechanisms considered in this paper could create a rather high cosmological fraction of the 
energy density of the relic gravitational waves at very high frequencies and gives an opportunity on 
investigating the high GWs spectrum by present and future detectors. Unfortunately at the lower part 
of the spectrum flew significantly drops down. Still the planned interferometers DECIGO/BBO could 
be sensitive to the predicted GWs. It is noteworthy that the mechanism of GWs generation suggested 
here kills or noticeably diminishes GWs from inflation by the redshift of the earlier generated GWs at 
the PBH (MD) stage. 
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